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ABSTRACT 

Growing Neutrino quintessence describes a form of dynamical dark energy that could 
explain why dark energy dominates the universe only in recent cosmological times. This sce- 
nario predicts the formation of large scale neutrino lumps which could allow for observational 
tests. We perform for the first time N-body simulations of the nonlinear growth of structures 
for cold dark matter and neutrino fluids in the context of Growing Neutrino cosmologies. 
Our analysis shows a pulsation - increase and subsequent decrease - of the neutrino density 
contrast. This could lead to interesting observational signatures, as an enhanced bulk flow 
in a situation where the dark matter density contrast only differs very mildly from the 
standard ACDM scenario. We also determine for the first time the statistical distribution 
of neutrino lumps as a function of mass at different redshifts. Such determination provides 
an essential ingredient for a realistic estimate of the observational signatures of Growing 
Neutrino cosmologies. Due to a breakdown of the non-relativistic Newtonian approximation 
our results are limited to redshifts z ^ 1. 
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1 INTRODUCTION 

The presently accepted cosmological model describes the evolution 
of the Universe by assuming the existence of two distinct forms of 
gravitating energy that largely dominate the total cosmic density, 
leaving to the particles described by the standard model of parti- 
cle physics only a small fraction of the energy budget. These two 
components are a family of non-relativistic massive particles (Cold 
Dark Matter, CDM hereafter) that source the gravitational insta- 
bility processes leading to the formation of galaxies and galaxy 
clusters, and a diffuse component with a negative equation of state 
(Dark Energy, DE hereafter) that accounts for the observed acceler- 
ated expansion ( ,Riess et al.|1998[|Perlmutter et al.|1999||Kowalski| 
|et al.pOOS] ) of the Universe during the last six billion years of its 
evolution. 

For the latter component, the standard model assumes a cos- 
mological constant A, i.e. a perfectly homogeneous form of energy 
whose density remains constant in time despite the expansion of the 
Universe. Although this simple choice reproduces with remarkable 
accuracy a large amount of observational data, it is faced by two 
fundamental theoretical problems concerning its fine-tuned value 
(known as the "fine-tuning problem") and the coincidence of its 
domination over cold dark matter (CDM) only at a relatively recent 
cosmological epochs (known as the "why now?" problem). In order 
to provide possible solutions to these issues, alternative scenarios 



for the cosmic DE have been proposed. In particular, a viable alter- 
native consists in identifying the DE with a scalar field cj) (the "cos- 
mon") dynamically evolving in a self-interaction potential ^(0), 
as proposed by |Wetteri ch ( 1988 ) and Ra tra & Peebles| ( |1988| ). 

For the former component, instead, several candidates can be 
considered (see jBertone et al.|200 5| for a review), including weakly 
interactive massive particles (WIMPS) as the neutralino ( Goldberg] 
1983 , Ellis et al.|1984^ ), arising in supersymmetric extensions of the 
standard model, or light scalars as the a xion (.Preskill et al.,1983( 
[Rosenberg & van Bibber|2000| ) or geon ( [Beyer et al.|2010| ). Mas- 
sive neutrinos, being an already known family of massive particles 
beyond the standard model would seem the most natural candidate 
for dark matter. They are, however, ruled out as the main contrib- 
utor to the dark matter density by observational constraints on the 
evolution of structure formation processes (see e.g. [Bond & Szalay[ 
|1983j). 

Nevertheless, even a small fraction of dark matter in the form 
of massive neutrinos could play a crucial role in the evolution 
of the Universe and in the late stages of large scale structure 
formation, in particular if a direct interaction between massive 
neutrinos and a DE scalar field is present, as will be discussed in 
the present paper. 

Cosmologies with an increasing neutrino mass stopping the 
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evoltuion of dark energy due to a coupling between neutrinos and 
a DE scalar field, known as the "Growing Neutrino" scenario, have 
been first proposed by |Amendola et al.||2008| ); [Wette rich (2007) 
with the aim to provide a viable explanation of the "why now?" 
problem. In such models, in fact, the reason behind the present 
quasi- static behavior of the DE density and equation of state re- 
sides in the direct cosmon-neutrino coupling. As a consequence of 
the interaction, the neutrino mass grows in time depending on the 
value of the cosmon. However, like in standard coupled DE models 
( |Wetterich| 1995[ |Amendola|2000| ), the interaction with relativistic 
species is suppressed for symmetry reasons. This implies that the 
coupling between massive neutrinos and the DE scalar field gets 
active only when neutrinos become non-relativistic. 

In the context of Growing Neutrino models, therefore, the neu- 
trino mass may be very small at decoupling, and grow into the eV 
range only in a rather recent cosmological epoch. A transition be- 
tween relativistic and non-relativistic neutrinos will then occur only 
at relatively recent redshifts, z ^ 5 — 10. At higher redshifts neutri- 
nos will free-stream and behave as relativistic particles. At z < 5, 
the time evolution of the neutrino mass changes the evolution of the 
cosmon itself via the Klein Gordon equation: 



b" + 2n(i)' + a 



dV 



(1) 



with pjy and = Wupiy the energy density and pressure of the 
neutrinos. 

As soon as this happens, the coupling to the neutrino fluid ef- 
fectively stops the further time evolution of the cosmon field, and 
from this time on DE effectively behaves as a cosmological con- 
stant. In other words, when neutrinos become non-relativistic, the 
scalar field potential ¥{(/)) acquires a correction due to the neu- 
trino coupling and the scalar field feels an effective potential with 
a slowly evolving minimum, where the DE scalar field stops its 
evolution. Since the effective potential Vef / differs from the true 
potential V^cj)) only when neutrinos become nearly pressureless. 
Growing Neutrino models account for the sudden recent dark en- 
ergy domination by relating it to a 'cosmological event', i.e. neutri- 
nos becoming massive enough to behave as a non-relativistic fluid. 

When this transition happens, the cosmon field moves away 
from its early scaling solution. Subsequently, it almost stops in the 
minimum of the effective potential. More precisely, the cosmon 
does not stop but oscillates around the minimum of the effective 
potential, leaving a characteristic oscillatory imprint also on the 
neutrino mass. In this paper we will describe for the first time in 
detail the effect that such an oscillating pattern leaves on structure 
formation. 

The cosmon-neutrino coupling /3 can be large (as compared 
to gravitational strength, see e.g. Fardon et al. 2004,^Bjaelde et al.| 
|2008[|Brookfield et al.|2006| ), such that even the small fraction of 
energy density in cosmic neutrinos can have important cosmologi- 
cal effects. The fifth force mediated by the cosmon is substantially 
stronger than the standard gravitational interaction. Its effect be- 
comes important only once the neutrino mass is sufficiently big for 
neutrinos to be non-relativistic. When this happens, neutrinos feel 
the presence of the fifth force and can collapse into stable bounded 
structures of the type described in Brouzakis et al. (2008 ); Berna^ 
Idini & Bertola mi (2009 ). 

It has been shown that neutrino lumps form at redshift Zn\ ~ 
1 — 2, when the neutrino fluctuations become nonlinear ( |Mota et aT] 
|2008[ |Wintergerst et al.|2010'| ). Furthermore, neutrino lumps signif- 
icantly affect only large scales ( |Mota et al.|2008| [Pettorino et al. 
|2010| ) with typical sizes in the range of a few 10 — 100 Mpc. 



A correct evaluation of the size and time dependence of the 
gravitational potential induced by neutrino lumps is essential (and 
yet difficult to achieve) for a comparison of Growing Neutrino 
models with current available data. Furthermore, in a Newtonian 
approximation the role of the gravitational potential is taken by 
/3 Scj), with Scj) the local fluctuation of the cosmon or quintessence 
field. While the gravitational potential remains much smaller than 
one even for highly nonlinear neutrino lumps, the cosmon induced 
potential 13 Scj) can reach values of order one, signalling a break- 
down of linear theory in the cosmon sector. The difficulty for the 
computation of the neutrino-induced gravitational potential comes 
from the fact that an extrapolation of the linear evolution violates 
rapidly even the most extreme bound which would result if all neu- 
trinos of the visible universe were concentrated in one spot ( Win-| 
[tergerst et al.|2010| ). In order to get a meaningful picture of Grow- 
ing Neutrino quintessence a statistical analysis of the neutrino lump 
abundances as a function of their mass is required. In order to 
achieve such a statistical understanding of neutrino lumps, it is nec- 
essary to rely on non-linear methods, such as N-body simulations, 
that can allow neutrinos to distribute, merge and diffuse with time 
and within different lumps. 

Both the size and the time evolution of the neutrino-induced 
gravitational potential are important features that need to be under- 
stood in order to draw a realistic picture of the Growing Neutrino 
scenario. It was shown in Pettorino et al. ( 2010 ) that the lump po- 
tential can have important effects on the propagation of photons and 
therefore influence the CMB -anisotropics via the integrated Sachs- 
Wolfe (ISW) effect, which may be observed by a modification of 
the CMB- spectrum at low angular multipoles. The neutrino gravita- 
tional potential also influences the correlation between matter and 
photon fluctuations for which the present observations yield val- 
ues larger than expected in the cosmological ACDM mod el (|Jain| 
&RaT ston 1999, Inou e & Silk|2006| [Rudnick et al.|2007l |Samal| 
et al. 12008 , Giannantonio et al. 2008 ). Moreover, a sudden increase 



in could reco ncile (Ayaita et al. 2009 ) the presently observed 
large bulk flows ( [Watkins et al.|2009| ) with observational bounds on 
the matter fluctuations at similar scales ('Reid e t al.|2010| ). The size 
of the gravitational potential also determines the counter-effect of 
neutrino lumps onto dark matter structures, possibly erasing bary- 
onic oscillations from the CDM power spectrum if the neutrino- 
potential exceeds substantially a value around 10~^ ( Brouzakis 
|etal.|2011| ). 

As already mentioned above, the difficulty in providing a re- 
alistic estimate of the neutrino gravitational potential resides 
in the fact that linear theory breaks down. A first estimate of 
based on extrapolations relying on a careful matching between 
linear equations for /c < 10~^ h/M pc and nonlinear resul ts for 
k>2 - 10"^ h/Mpc was attempted in Pettorino et al. |( |2010| ), pro- 
viding also a first attempt to address the backreaction of small neu- 
trino lumps onto larger-scale lumps. 

The evolution of isolated lumps was addressed within 
non-linear hydrodynamic equations in Wintergerst et al. (2010), 
showing how they mimic very large cold dark matter structures, 
with a typical local gravitational potential 10~^ for a lump size of 
10 Mpc. Spherical collapse within Growing Neutrino quintessence 
( [Winterge rst & Pettorino|20T0| first showed that the extrapolated 
linear density contrast at collapse manifests a characteristic 
oscillating pattern that reflects the evolution of the background 
neutrino mass rriu (0) as the cosmon moves around the minimum 
of the effective potential. Spherical collapse cannot, however, take 
into account effects due to velocity directionality, which, as shown 
by |Baldi et aT] ( |2010l ), can be very relevant in the presence of a 
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fifth force and can only be addressed in the context of N-body 
simulations. 

In this paper we present for the first time the results of N-body 
simulations for Growing Neutrino models. This allows us to fol- 
low not only the growth of neutrino perturbations in the nonlinear 
regime, but also, for the first time, to estimate the statistical abun- 
dance of neutrino lumps as a function of their final mass. N-body 
simulations further permit to follow the initial stages of the forma- 
tion of neutrino lumps with a detail never reached before: CDM 
structures first form, as in the standard scenario; then, we directly 
show how neutrinos start to concentrate along CDM filaments, the 
latter seeding the subsequent growth of neutrino lumps. 

Important differences need to be taken into account with re- 
spect to CDM simulations. First, it is necessary to switch between 
a phase in which neutrinos are relativistic and feel no coupling to 
a phase in which neutrino fluctuations grow under the effect of the 
fifth force; this transition happens around z ~ 4 for the Growing 
Neutrino model presented in |Mota et aL] p008| ) and investigated 
here. Second, the strength of the coupling can be order 10^ stronger 
than gravitational attraction, with a consequent fast growth of neu- 
trino lumps which will deserve particular care, both for its impact 
on the computational cost of the simulations and, physically, for 
its consequences on particle velocities. Indeed, the latter eventually 
limit the redshift down to which it is possible to rely on N-body 
simulations, valid as long as velocities are small enough compared 
to the speed of light. We discuss this issue in detail for the first time 
and show that the formation of neutrino lumps follows an oscillat- 
ing pattern, depending on the background oscillations and on the 
interplay of mass and velocity dependent terms. 

After a rapid overview of Growing Neutrino quintessence both 
at the background and at the linear levels in Sects. [2] and [3] we 
summarize what we already know about neutrino lumps in Sect. [3] 
We then describe the set of N-body simulations carried out for the 
present work and discuss in Sect. |4] the regime of validity of our 
analysis, ultimately limited by neutrino velocities growing up to 
relativistic values. In section [5] we present our results: we are able 
to follow neutrino particles and their velocities. Neutrino lumps, 
first seeded by CDM filamentary structures, evolve in an oscillating 
manner, becoming more and less concentrated as time goes by. We 
provide a physical intuitive explanation for this behavior in Sub- 
sec. |5^ and we compute the halo abundance of neutrino lumps as 
a function of their mass in Subsec. [54l In Subsec. l531 we consider 
effects on the CDM matter power spectrum and give an estimate of 
the power enhancement due to the presence of neutrino structures. 
Finally, in Sect. [6] we draw our conclusions. 



2 GROWING NEUTRINO COSMOLOGIES 

Growing Neutrino models involve a coupling between the DE 
scalar field (cosmon) and massive neutrinos. At the background 
level, the expansion of the universe obeys the Friedmann equation: 

7^^=(^]' = ^^^p. (2) 



where primes denote derivatives with respect to conformal time r. 
The sum is taken over all components a of the energy density in 
the universe, including CDM, DE, neutrinos, baryons and radia- 
tion. The time evolution of the energy density pa for each species 
involves the equation of state Wa = Pa/pa- We use units in which 
the reduced Planck mass Mpi/Stt is set to one. 



A crucial ingredient in this model is the dependence of the 
average neutrino mass on the cosmon field 0, as encoded in the 
dimensionless cosmon-neutrino coupling /3, 

dlnrriu 



(3) 



We consider here a constant /3 < such that for increasing the 
neutrino mass increases with time 



(4) 



where fhu is a constant. For a given cosmological model with a 
given time dependence of 0, one can determine the time depen- 
dence of the neutrino mass mu{t). For three degenerate neutrinos 
the present average value of the neutrino mass mu{to) can be re- 
lated to the energy fraction in neutrinos (h ^ 0.72) 

3mjy{to) 



(5) 



94eVh^ ' 

For the purpose of this paper we will consider the model de- 
scribed in Mota et al. (2008 ), with a constant coupling /3 = —52. 
This choice corresponds to a present neutrino mass of about 2.1 
eV. In general, /3 can be a function of 0, as proposed in |Wetterich| 
P007| ) within a particle physics model, modifying equation (|4]) but 
leading to similar effects. We refer to |Mota et aL] (2008 ) for fur- 
ther details on this model and recall here for convenience only its 
essential ingredients. 

The cosmon exchange mediates an additional attractive force 
between neutrinos of strength 2/3^. The case /3 ^ 1 corresponds to 
a strength comparable to gravity. The dynamics of the cosmon can 
be inferred from Eq. ([T]|, where we choose an exponential potential 
|Wetterich||T988t|Ratra & Peebles | ( [T988l ) |Ferreira & Joyce] ( [T998l l: 

y(0)(xe-"^ . (6) 

The constant a is one of the free parameters of our model and de- 
termines the slope of the potential and thus the DE fraction at early 
times. Current bounds constrain it to be of the order a ^ 10 or 
bigger ( [Doran et al.|2007| ) and we assume a = 10 in the present 
work. The naturalness of the exponential potential in the presence 
of quantum fluctuations and a cosmon-matter coupling has been 
discussed by |Wetterich| ( |2008| ). 

The homogeneous energy density and pressure of the scalar 
field are defined in the usual way as 



. (7) 



We choose the value of 7 = —/3/a such that we obtain the correct 
present dark energy fraction according to the relation { [Amendolal 
let al.|2008| ): 

' a 



(8) 



16eV ' ' 

Dark energy and growing neutrinos follow the coupled conserva- 
tion equations: 

p^ = -31-L{l + Wct,)pcf, + I3(l)'{l-3wu)pu , (9) 

p'^ = -3n{l + Wu)pu-P(|)'{l-3w^.)p^. . (10) 

The sum of the energy momentum tensors for neutrinos and DE is 
conserved, but not the separate parts. We neglect a possible cosmon 
coupling to CDM, so that pc = —3'Hpc- 

Given the potential ([6]), the evolution equations for the dif- 
ferent species can be numerically integrated, providing the back- 
ground evolution shown in Fig.[T]for a constant /3 = —52. The ini- 
tial pattern is a typical early dark energy model, since neutrinos are 
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l + z 

Figure 1. Energy densities of neutrinos (dashed, red), cold dark matter 
(solid, black), dark energy (dotted, green) and photons (long dashed, blue) 
are plotted vs redshift. For all plots we take a constant (3 = — 52, o; = 10 
and present neutrino mass = 2.11 eV. 



still relativistic and almost massless, with pi, = yOi//3, so that the 
coupling term in eqs.([TJ, (|9]), |T0| vanishes. DE is still subdominant 
and follows the attractor solution provided by the exponential po- 
tential (see |Wetterich|1995l|Amendola|2000l|Copeland et al.|1998[ 



Mota et al. 



2008 



for details): for 2; ^ 6 it tracks the dominant back- 
ground component with an early dark energy fraction = n/a^ 
and n = 3(4) for the matter (radiation) dominated era. As soon as 
neutrinos become nonrelativistic, the coupling term ^ I3p u in the 
evolution equation for the cosmon (Eq.[TJ starts to play a significant 
role. This cosmological "trigger event" kicks out of the attractor. 
As shown in Fig.[T] and change behavior for z < 6: the value 
of the cosmon stays 'almost' constant and the frozen cosmon po- 
tential mimics a cosmological constant. According to this model, 
at the present time neutrinos are still subdominant with respect to 
CDM, though in the future they will take the lead. For our choice 
of the parameters neutrino pressure terms may be safely neglected 
for redshifts z^r < 4; before that redshift neutrinos free stream as 
usual relativistic particles. 

Most importantly, the — zy coupled fluid manifests small 
decaying oscillations: the cosmon moves around the minimum of 
the effective potential and pu oscillates towards an almost constant 
value. The cosmon background oscillations enter also in the evo- 
lution of the neutrino mass miy{(j)) which in turn influences the in- 
duced gravitational potential. Consequently, the extrapolated linear 
density at collapse oscillates in redshift, as shown in |Wintergerst &| 
|Pettorino| ( [20T0l ). 

In this paper we will address the consequences of such a pat- 
tern in great detail, showing for the first time how these oscilla- 
tions directly affect structure formation. Before moving to the re- 
sults of N-body simulations, however, we recall in the next sec- 
tion our present knowledge about large scale structures in Growing 
Neutrino models and, in particular, on "neutrino lumps". 



3 NEUTRINO LUMPS 

Linear perturbations in Growing Neutrino quintessence have been 
widely described in |Mota et al.| ([2008) to which we refer for details. 
We recall here that, as it happens whenever a coupling of this type is 
active ( [Amendola|2004||Pettorino & Baccigalupi,2008, Baldi et aL| 
|2010| ), it is possible to reformulate the Euler equation in order to 
include the coupling contribution for each particle a: 

< + {n- 13(1)') vc. - V [$c. + = . (11) 

In cosmic time (dt = a dr) Eq. |TT| can also be rewritten in the 
form of an acceleration equation for particles at position r: 



-i/v„ 



, Ga,7Tla 



(12) 



The acceleration equation |T2| contains all the key ingredients that 
affect structure formation and it is therefore important to keep it in 
mind for a better intuitive understanding on the results presented in 
the next sections. These ingredients are: 

(i) a fifth force V + with an effective = Gn[1 + 
2/3^] ; 

(ii) a velocity dependent term |3(j)^/a which may act as dragging 
or friction according to the sign of 0, replacing effectively H 

(iii) a time-dependent mass for each particle a, evolving accord- 
ing to (|4]). 

We use here the cosmological background mass, while for a more 
accurate treatment rria would also be space dependent according to 
the local value of 0. The same holds for 0. A relativistic general- 
ization of eq.([TTJ is presented in the appendix but not used in the 
present paper. 

In a Newtonian approximation /3 6(1) exceeds the gravitational 
potential by a factor 2/3^ 5 x 10^ in the model under investi- 
gation. A local value /3 ^0 = 1 means that the local mass of the 
neutrinos is smaller by a factor as compared to the cosmolog- 
ical neutrino mass. As shown in Mota et ah] p008| ); |Wintergerst| 
|et al.| ( [20T0| ); |Pettorino et al.| ( [2010| one expects substantial modi- 
fications of linear theory once P 5(t) reaches values of order unity. 
Indeed non-relativistic neutrinos cluster under the effect of the fifth 
force, at scales estimated to be around a few 10 — 100 Mpc ( |Mota| 
|et al.|2008| ). At the end, they may form stable neutrino lumps of the 
type described by |Brouzakis et al.jp008| ); |Bernardini & Bertolami| 
(2009). 

When perturbations become non-linear, at a redshift z ^ 1 — 2, 
difficulties in evaluating the characteristic features of the lumps in- 
crease. Yet, it is crucial to estimate the size and time evolution of 
the gravitational potential of the neutrino lumps. The gravitational 
potential influences the low angular momenta of the CMB via the 
ISW ( |Pettorino et al.||2010| ) and determines the size of the cross- 
correlation between CMB and LSS. As a general rule of thumb, 
cosmologically avaraged potentials larger than ^ 10~^ at length 
scales of 100 Mpc or more give too strong an effect on the CMB. If 
substantially bigger than 10 the gravitational potential can also 
potentially erase baryonic acustic oscillations from the CDM power 
spectrum, thus constraining the range of values allowed for the cou- 
pling to match current available data dBrouzakis et al.||2bll| ). On 
the other hand, potentials around 10~ or somewhat smaller may 
cause an observational ISW effect, perhaps even with oscillatory 
structure. Furthermore, large neutrino lumps and a rapidly evolv- 
ing may correspond to large peculiar velocities, possibly recon- 
ciling (as suggest ed by Ayaita et al . 2009) the presently observed 
large bulk flows fWatkins et al.|2009| ) with observational bounds on 
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the matter fluctuations at similar scales ( |Reid et al.|2010| ). Growing 
Neutrino scenarios were also shown to be one of the first exam- 
ples of a realistic large backreaction mechanism which modifies 
the parameters of the homogeneous and isotropic background field 
equations by order one effects ( [Pettorino et al.|20ro[ |Nunes QiH] 



A first study of the gravitational potential for individual neu- 
trino lumps was performed in |Wintergerst et al.| p010| ) solving 
non-linear Navier-Stokes fluid equations for isolated lumps. This 
method allowed to follow the initial growth of non-linear fluctua- 
tions in physical space, by integrating numerically on a 3D grid 
nonlinear evolution equations, until virialization naturally occurs. 
Neutrino lumps were shown to form and mimic very large cold 
dark matter structures, with a typical local gravitational potential 

^ 10~^ for a lump size of ^ 10 Mpc, reaching larger values 
for larger lumps. 

In P ettorino et al.| ( [20T0l ) it was shown that a simple extrapo- 
lation of linear perturbations similar to [Franca et al.| ( [2009] ) to red- 
shift z = would incorrectly overcome even the extreme bound 
in which all available cosmic neutrinos collapse within one single 
lump; this limit value reaches at most ^ 10~^ for scales of 
about /c - 3 X 10"^ Mpc"\ In a more realistic scenario, in which 
neutrinos distribute within different lumps and merging is active, a 
substantially lower value of is expected. 

In the absence of accurate results for ^^{k), sl first estimate 
of the plausible overall size of and the ISW effect was per- 
formed in jPettorino et al.| ( |2010| ) by imposing reasonable bounds 
on the growth of fluctuations, with a careful matching between lin- 
ear and non-linear analysis. As no information on the statistical 
abundance of neutrino lumps was available, non-linear estimates 
had to rely on reasonable guesses on the neutrino statistical distri- 
bution in more lumps of one or two fixed typical sizes. Interest- 
ingly, these analysis had shown some hints of oscillating patterns 
possibly affecting the gravitational potential. The ISW-effect on the 
CMB power spectrum hinted to the presence of oscillatory features, 
reflecting oscillations in the cosmological neutrino energy density 
as visible in Fig. [T] One may speculate that this could be respon- 
sible for small "enhancements" and "depressions" in the observed 
CMB spectrum at low multipoles, though the size and details of the 
neutrino-induced ISW-effect depend on the particular model and 
require a more reliable estimate of the gravitational potential. 

In Wintergerst & Pettorino |p010| ) spherical collapse was also 
attempted, showing the presence of a clear oscillatory pattern in 
the extrapolated linear density at collapse. Non-linear power spec- 
tra of CDM and neutrino perturbations in Growing Neutrino mod- 
els were also shown in |Brouzakis et al.| { |2011| ), using the approach 
of IPietroml ( |2008| ), named Time Renormalization Group, or TRG. 
Such non-linear quantitative treatments cannot, however, account 
for a slowing down of the growth as compared to the linear approx- 
imation and only temptative reasonable bounds can be formulated: 
large uncertainties remain when exploring the region in which the 
neutrino sector becomes highly non-linear (z < 2). 

In the next sections we present for the first time results from N- 
body simulations of Growing Neutrino cosmologies, following the 
evolution of neutrino lumps and their merging through cosmic his- 
tory until z ^ 1. We show that the evolution of neutrino structures 
is initially seeded by CDM filaments, that form first as in a stan- 
dard cosmological scenario. Then, the growth of neutrino lumps 
indeed proceeds in a clear oscillatory fashion, more pronounced 
than ever expected in previous works. As we will see in the fol- 
lowing, this is also due to the fact that N-body simulations finally 
allow us to take into account the effects due to the directionality 



of velocity-dependent terms like the one in eq.|T2|, not included 
within a spherical collapse picture: switching from a 'friction' -type 
to 'dragging' -type term is determinant to enhance the oscillating 
evolution of neutrino lumps. Furthermore, we are able for the first 
time to evaluate the statistical abundance of neutrino lumps as a 
function of their mass. 



4 N-BODY SIMULATIONS OF GROWING NEUTRINO 
COSMOLOGIES 

In order to study the nonlinear evolution of structures in the context 
of Growing Neutrino cosmologies, we run a series of N-body sim- 
ulations that take into account the interaction of neutrinos with the 
cosmon and all the related effects on the growth of density perturba- 
tions that have been discussed in the previous Sections. To this end, 
we have made use of the modified version by "Ba ldi et ar](|2010| ) 
of the cosmological N-body code GADGET-2 ( Sprin gel|2005| ), that 
was specifically designed to include the effects of interacting DE 
models. 

The only modification that has to be implemented in the code 
with respect to the algorithm described in |Baldi et al.| ( |2010] | con- 
cerns the transition of the neutrino particles from the relativistic 
to the non-relativistic regime. Due to the fast growth of the neu- 
trino mass, this transition can be roughly approximated with a step 
function for the neutrino equation of state, sharply changing from 

= 1/3 to = at the transition redshift Znv In this study 
we have assumed Znr = 4 based on the results of the background 
evolution of the specific Growing Neutrino model under investiga- 
tion described in Sec.|2] According to this approximation, neutrino 
particles in the simulation box are free-streaming at all scales for 
z > Znr, when their mass is negligible with respect to the CDM par- 
ticle mass, while they behave as cold dark matter for z < Znr when 
their growing mass rapidly makes them non-relativistic. To imple- 
ment the transition in our numerical machinery, we have artificially 
switched off the gravitational interaction for neutrino particles be- 
fore Znr, and set their mass to zero. Therefore, the neutrino particles 
in the simulation box will neither experience any acceleration nor 
exert any force on the CDM particles during their relativistic regime 
z > Znr, and will therefore move on straight trajectories with their 
own initial velocities. The initial velocities are set in the simulations 
initial conditions by randomly sampling a Fermi-Dirac distribution 
with a temperature corresponding to the average neutrino tempera- 
ture at Znr according to our model of growing neutrino mass. This 
implementation ensures that the neutrino density field will remain 
completely homogeneous in the simulation box until Znr, and that 
the velocity distribution of the neutrino particles at Znv corresponds 
to the correct Fermi-Dirac velocity distribution for a neutrino with 
mass mu^Znv) ~ 0.02 eV. 

After the transition redshift Znv, neutrino masses are set to 
their correct value, according to the neutrino mass evolution of the 
specific model under investigation, and the gravitational interaction 
as well as the scalar fifth-force are switched on. Therefore, right 
after Znv neutrino particles will start to experience gravitational 
accelerations sourced by the surrounding CDM density field, that 
by Znv will have already developed significant inhomogeneities. 
Furthermore, neutrino particles will also start exerting gravitational 
forces on other CDM and neutrino particles, due to their rapidly 
growing mass. It is however the fifth-force acting among neutrino 
particles that soon dominates over gravity, driving the large scale 
evolution of neutrino structures at 2; < Znr- 
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Model 


Box Size 
(Mpc / h) 


Neutrino 
Particles 


CDM 
Particles 


Neutrino part. 

mass at 2; = 
Uvi0/n; 


CDM part. 

mass 
Uvi0/n; 


Neutrino grav. 
softening 


CDM grav. 
softening 


ACDM-a 


40 


643 


1283 




2.0 X 10^ 




8.0 


Cnu-a 


40 


643 


1283 


8.7 X 10^ 


2.0 X 10^ 


15.0 


8.0 


ACDM-b 


120 


643 


1283 




5.4 X 10^0 




24.0 


Cnu-b 


120 


643 


1283 


2.3 X IQii 


5.4 X 10^0 


45.0 


24.0 


ACDM-c 


320 


1283 


2563 




1.3 X 10^1 




16.0 


Cnu-c 


320 


1283 


2563 


5.6 X IQii 


1.3 X 10^1 


30.0 


16.0 



Table 1. The different simulations considered in the present work for Growing Neutrino cosmologies and for the standard ACDM model. 



With this implementation we have run a series of N-body sim- 
ulations of structure formation for the specific Growing Neutrino 
model described in Sec. [2] These consist in cosmological boxes of 
different sizes containing both neutrino and CDM particles. The 
initial conditions for the CDM component have been generated by 
setting up a random-phase realization of the Eisenstein & Hu power 
spectrum (Eisenste in & Hu|1997| ) according to the Zel'dovich ap- 
proximation (Zel'dovich 1970 ). The normalization amplitude of the 
power spectrum is adjusted to a value of as = 0.769. The displace- 
ments of the particles are then rescaled to the starting redshift of 
the simulation Zi — 60 with the linear growth factor which is 
assumed to be the ACDM one. Due to the impact of the growing 
neutrino component on the growth of structures, however, the real 
CDM growth factor might slightly differ from the standard ACDM 
one, thereby driving to an actual value of as different from the one 
used for the normalization of the initial conditions. Since the neu- 
trino component is assumed to be completely homogeneous before 
Znr, no displacement is applied to neutrino particles. The procedure 
adopted to generate the initial velocities for the neutrino component 
has already been described above. 

A major numerical challenge for the study of Growing Neu- 
trino cosmologies by means of N-body simulations consists in the 
significant increase of the overall computational time compared 
to standard CDM simulations with the same number of particles, 
which can be as large as a factor of 10^. This large increase of 
the required computational time is mainly due to the timestepping 
criterion, which relates the individual timestep of a given parti- 
cle to the acceleration experienced by the particle in the previous 
timestep, according to the equation: 

At, oc Ve/|aoid| , (13) 

where e is the gravitational softening and add is the gravitational 
acceleration of the particle. With the above timestepping criterion, 
as soon as the neutrino fifth-force becomes active right after Znr, 
the acceleration experienced by neutrino particles exceeds by a 
factor 2/3^ the standard gravitational acceleration, and the corre- 
sponding timestep At is therefore reduced by a factor 1 / 1/2^, 
that for the model under investigation is of the order of ^ 80. This 
in turn determines a similar increase in the total computational 
time of the simulations. The situation is well illustrated in Fig. [2] 
where the cumulative wallclock time (in seconds) for two of our 
simulations is plotted as a function of the cosmological scale 
factor a. The solid red line represents the simulation for coupled 
neutrinos, while the black dashed line is for the standard ACDM 
cosmology. As one can clearly see, while before the transition 
redshift Zr^x the two simulations evolve with identical wallclock 
times, after the neutrino fifth-force becomes active at 2; = Znr the 
coupled neutrino simulation starts deviating from the ACDM one, 
with a rapidly increasing computational cost. In order to alleviate 
this problem, since the increase of the total computational time for 
coupled neutrino simulations arises from the reduced timestep of 



E 





1 1 1 1 1 1 1 1 1 1 1 1 1 1 _ 
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Znr ACDM 
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a 

Figure 2. The cumulative wallclock time (in seconds) as a function of the 
cosmological scale factor a for a standard ACDM simulation (black dashed 
line) and a Growing Neutrino simulation with the same box size and particle 
number. The vertical line indicates the transition redshift Znr at which the 
coupling between DE and massive neutrinos becomes active. 



neutrino particles, and since neutrino structures are expected to 
form only at relatively large scales (as predicted by e.g. |Mota et al.| 
|2008| ) thereby requiring a relatively low spatial resolution, we have 
limited the number of neutrino particles in each simulation to be 
8 times smaller than the corresponding number of CDM particles. 
This contributes to reduce the increase of the total computational 
cost of the simulations, resulting in a relative overhead of a factor 
of 25 for our coupled neutrino runs as compared to ACDM (see 
again Fig. [2]), significantly smaller than the overhead of a factor of 
80 estimated for an equal number of CDM and neutrino particles. 

The details of the simulations considered in the present work 
are described in Table [T] 



5 RESULTS 

We now describe the results obtained from the analysis of the dif- 
ferent simulations described in Table[T] First of all, we focus our at- 
tention on the range of validity of the N-body approach for Growing 
Neutrino cosmologies, then we study the evolution and the statisti- 
cal properties of neutrino structures at different scales, and finally 
we investigate the backreaction effects of the growth of neutrino 
structures on the CDM distribution. 
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Figure 3. Left panel: The evolution of the neutrino mass rriu /mu{tQ) (solid red line) and of the average neutrino peculiar velocity in units of the speed of 
light < (Tu > (black line, squares) as extracted from the Cnu-b simulation described in Table[2 as a function of redshift. Right panel: the fraction of neutrino 
particles in the Cnu-b simulation of Table [T] exceeding a fraction X of the speed of light as a function of redshift. The different curves are for X = 0.1 (soHd 
black), X = 0.2 (dotted red) and X = 0.3 (dashed green). 



5.1 Evolution of neutrino velocities 

As a first step in our analysis we want to test the validity of the 
N-body approach to study the nonlinear evolution of structure for- 
mation in the context of Growing Neutrino cosmologies. N-body 
codes are designed to follow the evolution of density perturba- 
tions under gravitational instability processes by solving the Pois- 
son equation in the newtonian limit. This means that the scales 
of interest should be small compared to the cosmological horizon 
and that particle velocities should be small compared to the speed 
of light. A third issue concerns the field dependence of the neu- 
trino masses in our model. The cosmon field inside neutrino lumps 
may differ from its cosmological value, and as a consequence the 
neutrino mass may depend on position. This "backreaction" effect 
could be substantial for very large neutrino lumps (Pet torino et al.| 
|2010[[Nunes et al.|2011| ). Its computation would require to follow 
the scalar field distribution during the simulation, which goes be- 
yond the scope of the present implementations of N-body codes. 
(In our approach we use the cosmological background value of the 
cosmon for the determination of the neutrino mass, such that rriu 
depends on time but not on position). 

For the first issue, |Chisari & ZaldarriagaT2011) have recently 
shown that even for scales comparable or larger than the cosmolog- 
ical horizon the standard N-body algorithms provide a correct com- 
putation of the gravitational potential. With our simulations we can 
now address the second issue related to particle velocities. Particle 
velocities constitute a potential problem for N-body codes since in 
the newtonian limit particles are allowed to accelerate to arbitrarily 
large values, possibly even exceeding the speed of light. While this 
is not a problem for standard cosmological simulations, where pe- 
culiar velocities are always small with respect to the speed of light 
due to the smallness of the newtonian gravitational constant Gn , 
it might turn out to be a problem for non-standard scenarios - as 
e.g. the Growing Neutrino model under investigation in this work - 
where the attraction is governed by an effective strength GeS that 
in some situations significantly exceeds the value of Gn • 

For the specific case of the model considered in the present 
study, while the gravitational dynamics of CDM particles is al- 
ways governed by the newtonian gravitational constant Cn, the dy- 
namical evolution of neutrino particles after the transition redshift 



Znv = 4 is driven by an effective coupling Geff = G^{1 + 2/3^), 
which is more than 5 x 10^ times larger than Cn- It is therefore 
natural to ask whether this large effective gravitational constant 
leads to an acceleration of neutrino particles to exceedingly high 
velocities, and to which extent an N-body treatment of this sys- 
tem can be considered correct. In order to investigate this problem, 
we compute the evolution of the average neutrino peculiar veloc- 
ities by taking the mean of the neutrino particles velocities in our 
simulation boxes at different timesteps. The result of this procedure 
for one of our simulations (the Cnu-b run of Table [T]) is shown in 
the left panel of Fig. [3] where the solid black line and the black 
squares represent the average neutrino peculiar velocity extracted 
from the snapshots of the simulation in units of the speed of light 
< (Jzy >, where = Vu/c. As expected, we find that neutrino 
velocities rapidly grow after the transition redshift Znv reaching al- 
ready at 2; ^ 1 a non-negligible fraction (20 — 30%) of the speed 
of light, and exceed the value of c at 2; ^ 0.5. The final value of 
the average neutrino velocity in the simulation at z = is almost 
two times larger than the speed of light, having reached clearly a 
regime where N-body simulations are no longer reliable. 

A comparison with the non-linear analysis done in I Winterg^ 
[erst et al.] ([2010 ) solicited by the present study has confirmed the 
same trend of neutrino velocities shown in the left panel of Fig. |3] 
(N. Wintergerst, /?nva^^ communication). 

It is therefore clear that the N-body approach is not sufficient 
to follow correctly the dynamical evolution of neutrino particles 
up to the present time, and that relativistic corrections would be 
required in order to implement in the treatment the concept of a 
speed limit. This however would go beyond the newtonian limit 
in which N-body codes are implemented, and would require the 
development of completely new algorithms. In the present work, 
we do not want to undertake such a challenging enterprise, but 
we will rather try to assess the range of validity of the standard 
N-body approach and will restrict our analysis to this limited 
range. In particular, as a result of our analysis shown in Fig. [3] we 
can assume our simulations to be no longer reliable at redshifts 
below z ^ 1. We will therefore discard all the simulations outputs 
at 2; < 1 and will limit our analysis to the early stages of the 
nonlinear evolution of neutrino large scale structures. 
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The other interesting feature emerging from the analysis 
shown in the left panel of Fig. [3] is that neutrino velocities do not 
grow monotonically in time, but rather exhibit an oscillatory be- 
havior corresponding to a series of subsequent accelerations and 
decelerations of neutrino particles. This peculiar evolution is not a 
numerical artifact, but is directly related to the evolution of neu- 
trino masses. As discussed in Sec. [2] when neutrinos become non- 
relativistic the coupling with the cosmon gets active and the effec- 
tive potential Kff for the scalar field acquires a minimum where 
the scalar field stops, thereby mimicking a cosmological constant. 
However, due to its initial kinetic energy, before stopping the scalar 
field oscillates around the minimum of Feff, and these oscillations 
determine a corresponding oscillation of the neutrino mass accord- 
ing to Eq. |4] 

In the left panel of Fig. [3] the solid red line represents the 
evolution of the neutrino mass mu/mi^(to) ^ exp(— /3(/)/Mpi). 
As it clearly appears from the plot, the oscillations of the average 
neutrino velocity < au > are correlated with the oscillations in 
the neutrino mass. This result can be easily interpreted by having a 
look at the neutrino acceleration equation, Eq.[T2l 

Furthermore, we may also infer from Fig. |3] that even for a 
large present average neutrino mass (to) = 2.1 eV, the neutrino 
mass at z = 1 is rriu ^ 0.3 — 0.4 eV, with even smaller values 
for 2; > 1. It is yet to be investigated whether this value of mu{to) 
remains compatible with restrictions from structure formation { |Ko-| 
Imatsu etal.|2011| ). 

The right panel of Fig.|3]shows - for the same simulation Cnu- 
b - the fraction of neutrino particles whose velocity exceeds a given 
fraction X of the speed of light. As one can see from the plot, 
almost 90% of particles (solid black line) have velocities larger than 
10% of the speed of fight at z 1, while 20 - 35% of particles 
(dashed green line) exceed 30% of the speed of light in the same 
redshift range. 

5.2 Large scale neutrino structures 

Having in mind the limitations of the N-body approach for the 
study of coupled neutrino cosmologies discussed in the previous 
Section, we now start to investigate the growth of neutrino and 
CDM structures in our simulations. We begin by visualizing the 
evolution of the neutrino and CDM particle distributions in all 
the simulations described in Table [T] for three different redshifts 
above our limiting redshift z — 1 and below the transition redshift 
Znv = 4. Fig.|4] Fig. [5] Fig. [6| show a random subsample of CDM 
(grey points, top rows) and neutrino (colored points, bottom rows) 
particles at redshifts 2; = 1.8 (left columns), 2; = 1.3 (central 
columns), and z = 1.0 (right columns) for our coupled neutrino 
simulations Cnu-a (Fig. |4]), Cnu-b (Fig. |5]) and Cnu-c (Fig. [6|. 
In these figures, neutrino particles have been colored according 
to their velocities, with blue points representing "slow" particles 
(cTi^ < 0.1), orange points representing "intermediate velocity" 
particles (0.1 < < 0.3) and red points representing "fast" 
particles (a^ > 0.3). 

By having a look at the Cnu-a simulation (Fig. |4]) one can 
clearly see how CDM structures (first row) evolve on scales of the 
order of few tens of Mpc. The main filamentary structure is already 
present at 2; = 1.8 and evolves by growing overdensities and en- 
larging voids up to 2; = 1. It is very interesting to see what happens 
to the neutrino distribution (second row) at the same scales and 
in the same redshift range. At 2; = 1.8 neutrinos have a mass of 
^ 0.2 eV and had sufficient time after the transition redshift Znv to 



fall into the potential wells of the already developed CDM struc- 
tures. In fact, as one can clearly see in the figure, the distribution of 
neutrino particles at 2; = 1.8 roughly follows the underlying CDM 
distribution, with the main concentrations of neutrinos located in 
the same positions as the underlying CDM halos and filaments. At 
this epoch, however, neutrinos did not have enough time to acceler- 
ate to very high velocities, and in fact the figure shows that basically 
none of the neutrino particles has a velocity larger than 0.1. 

The situation is already significantly different at 2; = 1.3, 
where the first isolated neutrino lumps have formed, and where the 
distribution of neutrino particles does not trace anymore the un- 
derlying CDM filamentary structure, but rather concentrates in a 
few very large and isolated halos. At this time, some particles with 
higher velocities are clearly visible in the central regions of the 
largest neutrino lumps, while there is almost no diffuse neutrino 
component left over outside the main bound structures: practically 
all neutrinos are residing in bound halos at this redshift. 

The situation at 2; = 1 therefore appears at a first glance quite 
surprising: while the number of large neutrino halos in the box has 
been reduced to a few by the merging of a large number of smaller 
neutrino lumps, a significant diffuse component of free neutrinos, 
that was completely absent at 2; = 1.3, surprisingly appears again. 
We will find this phenomenon also in the larger scale simulations 
Cnu-b and Cnu-c, and we will come back later on its interpretation. 
The neutrino distribution at 2; = 1 in this small scale simulation 
shows the presence of a few neutrino lumps with a scale of 5 — 10 
Mpc/h in a volume of 6.4 x 10^(Mpc/h)^, and a dominant fraction 
of "slow" neutrinos, in agreement with Fig.|3] 

If we now consider the Cnu-b simulation (Fig. [5]) we can see 
how the CDM structures are significantly less evolved (as expected) 
on these larger scales as compared to the Cnu-a simulation. How- 
ever, the large scale filaments and the main CDM structures are vis- 
ible already at 2; = 1.8, and become progressively better defined as 
time goes by until 2; = 1. Also in this case, however, the neutrino 
distribution at 2; = 1.8 seems to follow the main features of the 
CDM structures, with neutrino particles distributed in very large 
scale lumps and filaments, and with still low velocities as basically 
all the neutrinos appear as "slow" particles at this stage. 

At 2; = 1.3 the situation is then completely different, with 
a large number of neutrino structures at very different scales, and 
with the presence of "slow" neutrinos as well as "intermediate" 
and "fast" neutrinos. As already found for the Cnu-a simulation, 
also in this case most of the neutrinos are residing in lumps, and 
there is basically no diffuse component of free neutrinos at this 
stage, i.e. the space between neutrino lumps is practically empty of 
neutrinos. 

Once again, the situation changes at 2; = 1, where the number 
of neutrino structures has significantly reduced, and one can 
identify neutrino halos on scales ranging from a few Mpc/h up 
to 20 — 30 Mpc/h. Also in this case, besides a few lumps that 
are still surrounded by "slow" neutrinos and a majority of lumps 
made by "intermediate" and "fast" neutrinos, a significant diffuse 
component of free neutrinos appears again, after being absent at 
z = 1.3. 

We can now move to consider the largest scale simulation at 
our disposal, the simulation Cnu-c (Fig. [6|. At these very large 
scales (320 Mpc/h aside), the CDM distribution (upper row) ap- 
pears very homogeneous (as expected) and no significant structures 
can be easily identified even at 2; = at our resolution. Neverthe- 
less, the CDM gravitational potential is sufficiently evolved, even 
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Figure 4. The distribution of CDM (top row) and neutrino (bottom row) particles in the Cnu-a simulation at z = 1.8 (left) z = 1.3 (middle) and z = 1 
(right). The neutrino particles are colored according to their velocity (Jiy, as explained in the legend. 



at these large scales, to source the growth of neutrino overdensi- 
ties right after the transition redshift Znr such that at z = 1.8 the 
neutrino distribution already shows a significant level of structure 
that is not visible in the underlying CDM distribution. This shows 
very clearly how the neutrino distribution evolves in the context 
of Growing Neutrino scenarios: at Znr, the neutrino component is 
completely homogeneous and suddenly starts to experience accel- 
erations sourced by the gravitational potential of the underlying 
CDM distribution; as soon as neutrinos, as a consequence of these 
accelerations, start to move towards the minima of the CDM gravi- 
tational potential and develop their own inhomogeneities, the fifth- 
force acting between neutrino particles and mediated by the scalar 
field 0, that is 5 x 10^ times stronger than gravity, starts driving the 
evolution of neutrino overdensities that can therefore grow much 
faster than the underlying CDM density field, and develop signif- 
icant structures at scales where CDM appears to be almost com- 
pletely homogeneous. The neutrino distribution at 2; = 1.8 in Fig.[6] 
therefore shows very clearly how the CDM density field acts as a 
"seed" for the growth of neutrino lumps. 

At this stage, also on these large scales, the neutrino com- 
ponent is almost completely made of "slow" particles, while at 
z = 1.3 a large number of isolated neutrino lumps on a wide 
range of scales have already formed and are mainly composed of 
"intermediate" and "fast" neutrinos. Once more, there is basically 
no diffuse component of free neutrinos at 2; = 1.3 and all neutrinos 
are residing in lumps, but the situation changes at z = 1, where 
a large number of neutrino lumps is still visible but where a 
significant diffuse component of free neutrinos appears again. In 
this case, the large majority of neutrino particles at 2; = 1 is "fast" 



while very few "slow" neutrinos are visible in the simulation box. 

The analysis of the evolution of the neutrino distribution in 
these three simulations of different scales that was described above, 
represented in Fig.|4] Fig.[5]and Fig. [6] clearly identifies two inter- 
esting and distinct phenomena. 

First, if we look at the neutrino distribution at z = 1 in the 
three simulations, we can notice a clear difference: the neutrino 
component is mainly composed of "slow" neutrinos in the Cnu- 
a simulation (40 Mpc/h aside), of "intermediate" neutrinos in the 
Cnu-b simulation (120 Mpc/h aside), and of "fast" neutrinos in the 
Cnu-c (320 Mpc/h aside). This indicates that the inclusion of larger 
scale modes in the CDM density field determines a larger integrated 
acceleration for the neutrino particles, since neutrinos have more 
time to accelerate in the large scale linear CDM potential than in 
the small scale nonlinear one. 

Besides possible questions of numerical convergence, the se- 
quence of different scales for the boxes indicates an issue for the 
scaling of the simulations. While we can find the structures of the 
small scale simulations also in the simulations at larger box size, 
the typical velocities for structures of a given size are not the same. 
In Fig|4]the largest structures have typical size of a few Mpc/h, 
with nonrelativistic velocities of the "particles". Structures of a 
similar size in Fig |5|6| have typically higher velocities. The rea- 
son seems simple: the lumps of say 5 Mpc/h size move as whole 
with a high peculiar velocity in the simulations with large box size, 
while they feel no "outside attraction" in the box of a compara- 
ble size and do therefore not develop a collective velocity. In the 
three simulations no saturation of this phenomenon is seen at the 
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Figure 5. The same as for Fig.|4]but for the Cnu-b simulation of Table[T] 
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largest box size. Investigating still larger box sizes is presumably 
not of much help since the peculiar velocities become so large that 
a nonrelativistic approximation seems no longer meaningful. It is 
not known how the collective velocities and scattering of smaller 
size structures influence precisely the neutrino distributions inside 
these structures. The conclusion drawn from our simulations should 
therefore only be considered as semi-quantitative. Nevertheless, a 
few interesting lessons can be drawn. 

5.3 Oscillating structure formation 

We have found in all the three simulations, independently on the 
size of the simulation box, that the neutrino distribution has a sig- 
nificant diffuse component at z = 1.8, then shows basically no 
diffuse component at z = 1.3, when all neutrinos are residing in 
lumps, and finally a significant diffuse component appears again 
at 2; = 1. This means that a large number of neutrino particles 
that were residing in gravitationally bound lumps at z = 1.3 are 
"released" and they appear no longer bound at 2; = 1. This ef- 
fect seems surprising since in standard gravitational instability pro- 
cesses for collisionless dynamics, bound objects remain bound and 
the diffuse component of particles continuously decreases. In other 
words, in the standard picture of gravitational instability density 
fluctuations can only grow in time. In this case, instead, we are 
witnessing an oscillatory behavior of the gravitational instability 
of neutrino overdensities that first collapse into a collection of iso- 
lated bound objects, separated by large empty regions, and then 
decay again releasing a significant fraction of their content into a 
new homogeneous diffuse neutrino field. 

This behavior can be understood by considering again the evo- 
lution of the neutrino mass and the corresponding oscillatory be- 
havior of the average neutrino velocity discussed in Sec. |5.1| In 
particular, the simultaneous time evolution of the neutrino mass and 
of the friction term appearing in the neutrino acceleration equation 
(Eq. [12]) can account for the results of our simulations. To better 
understand the situation, one can have a look at Fig. [7] where the 
average friction term < > /HMpi as extracted from our 
numerical simulations (top panel) and the neutrino mass evolution 
(bottom panel) are plotted as a function of redshift in the redshift 
interval 0.95 — 2. The filled red squares indicate the values of the 
two quantities in correspondence of the redshifts shown in Fig. [4] 
Fig.|5]andFig.[6] 

As a consequence of the oscillations of the cosmon around 
the minimum of its effective potential, the friction term I3(j)v also 
oscillates and changes sign during the evolution of the universe. 
This implies that the friction term actually acts as a drag term, 
i.e. accelerates particles in the direction of their motion, whenever 
positive, while it behaves as a real friction, i.e. decelerates parti- 
cles, when it is negative. The vertical dotted lines in Fig.[7]indicate 
the locations where the friction term changes sign. These points 
of sign inversion clearly correspond to maxima and minima of the 
neutrino mass evolution: whenever the neutrino mass grows in time 
(rhu > 0), the friction term is negative and acts as a real friction on 
neutrino particles, decelerating their motion, while when the neu- 
trino mass decreases in time (rhu < 0), the friction term becomes 
positive, acting as a drag term and accelerating particles. By hav- 
ing in mind the contextual evolution of the neutrino mass and of 
the friction term, it is possible to account for the results shown in 
Fig.|4lFig.|5]andFig.|t] 

Neutrinos evolve from a highly homogeneous distribution 
at 2; = Znr and start falling into the CDM potential wells as 
soon as they become non-relativistic. Therefore, the first neutrino 
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Figure 7. The evolution of the average specific friction term ^j) < > 
I HMpi as extracted from the Cnu-b simulation of Table[T](top panel), and 
the evolution of the neutrino mass in eV (bottom panel) as a function of 
redshift, in the redshift range 0.9 — 2. The red squares indicate the values 
of the two quantities at the redshifts represented in the different columns of 
Figs.[4]|5] and|6] i.e. 2; = 1.8, 2; = 1.3 and z = 1, respectively. 



structures follow the main filamentary pattern of the already 
evolved CDM distribution, as one can see at z = 1.8. However, 
as soon as some level of inhomogeneity has developed in the 
neutrino distribution, the scalar fifth-force acting only between 
neutrino particles overcomes standard gravity and becomes the 
main driver of the neutrino evolution. At this stage, i.e. between 
z = 1.8 and z = 1.3 the neutrino mass is constantly increasing, 
making the source term for the neutrino acceleration in Eq.[T2]also 
increasing. In other words, the neutrino potential wells become 
deeper and accelerate the instability process. At the same time, 
the friction term is negative, thereby acting as a real friction, 
i.e. decelerating particles. This reduces the neutrino kinetic energy 
thereby contributing to trap neutrino particles into gravitationally 
bound structures. This is what we can see in our snapshots at 
z — 1.3. After this time, however, the neutrino mass starts 
decreasing and consequently the neutrino potential wells start 
decaying, increasing the total energy of bound neutrino lumps. 
Contextually, the friction term has become positive and acts as 
a dragging force, accelerating particles in the direction of their 
motion, and therefore increasing the neutrino kinetic energy. These 
two phenomena contribute to increase the total energy of a large 
number of neutrino particles above zero, and to unbind them from 
the previously collapsed lumps, releasing a significant fraction of 
free particles that fill the space between the surviving neutrino 
halos. This is the situation that we observe in our snapshots at 
z = 1, where a significant diffuse neutrino component is present 
again in our simulation boxes. 

In conclusion, we have shown here for the first time that the 
oscillations of the scalar field around the minimum of its effective 
potential in Growing Neutrino cosmologies have a strong impact 
on the growth of linear and nonlinear neutrino structures, and de- 
termine a related oscillatory behavior of the evolution of neutrino 
density fluctuations. In other words, neutrino structures that form 
in the context of Growing Neutrino models do not continuously 
grow, but rather exhibit a sort of "pulsation" due to the oscillations 
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Figure 8. The cumulative neutrino halo mass function for the full sample of 
neutrino halos obtained by combining the separate neutrino halo catalogs 
of the three simulations Cnu-a, Cnu-b, and Cnu-c. Different curves show 
the neutrino halo mass function at different redshifts from z = 1.7 (red) to 
z = 1 (orange). The oscillation of neutrino structures is clearly visible in 
the inversion of trend for the neutrino halo mass function before and after 
z = 1.4. 

of the DE scalar field that drives the growth of the neutrino mass. 
As a consequence of this "pulsation", also the gravitational po- 
tential associated to neutrino lumps and large scale structures will 
not evolve monotonically, but rather follow an oscillatory behavior, 
with the time derivative i>zy of the gravitational potential induced 
by neutrino overdensities changing sign along with the oscillations 
of the cosmon cj). This non-monotonic evolution of the gravitational 
potential time derivative i>zy could significantly reduce the impact 
of Growing Neutrino models on the amplitude of the ISW effect 
for the low-multipole CMB angular power spectrum with respect 
to previous linear and nonlinear estimates (as e.g. |Pettorino et al.| 
|2010| ). A quantification of the ISW effect associated to the neutrino 
halo "pulsation" will therefore be necessary in order to test the vi- 
ability of the model and could provide a direct observational way 
to constrain the neutrino coupling. A non-monotonic evolution of 
the neutrino-induced gravitational potential could also alleviate the 
tension between observed bulk velocities and the ISW-correlations 
between temperature fluctuations and observed structures ( [Ayaita| 
letal.,2009^ 

5.4 Neutrino halo mass function 

As a next step in our analysis, we aim to quantify the number den- 
sity of neutrino lumps as a function of their total mass and inves- 
tigate how this number density evolves in time. We are therefore 
seeking to extract from our numerical simulations a "neutrino halo 
mass function", analogous to the standard CDM halo mass func- 
tion. To this end, we first need to build a catalog of neutrino struc- 
tures formed in our simulated cosmologies. We therefore identify 
neutrino groups with a Friends-of-Friends algorithm with linking 
length A = 0.2 X J where J is the mean interparticle spacing, and 
we consider only groups containing at least 200 particles. For each 
of these groups we compute the total group mass MfoF and the 
position of the center of mass. 



We compute our catalogs at different redshifts for the three 
Growing Neutrino simulations Cnu-a, Cnu-b, and Cnu-c, and we 
then combine them into a single neutrino halo mass function in or- 
der to maximally extend its mass range. To do so, we build a single 
catalog of all the halos identified in the three different simulations 
and we associate to each halo the corresponding number density 
computed from the mass function extracted from its original simu- 
lation. This new halo catalog is then ordered by halo mass, and the 
halos are binned in five logarithmic mass bins covering the whole 
mass range of the catalog. The number density in each mass bin is 
then computed by averaging the original number density of all the 
bin members. The result of this analysis is shown in Fig. [8] where 
the full neutrino halo mass function extracted from the combina- 
tion of the three simulations Cnu-a, Cnu-b, and Cnu-c is plotted for 
the five mass bins of our full catalog and for five different redshifts. 
This result represents the first estimate ever made for the abundance 
of neutrino lumps in a realistic realization of the nonlinear evolu- 
tion of a Growing Neutrino cosmology. 

Once more, as discovered by the visual inspection of the neu- 
trino distribution as a function of redshift discussed in Sec. |5.2| 
we can notice that the abundance of neutrino lumps does not grow 
monotonically in time. On the contrary, as Fig.|8]clearly shows, the 
number density of neutrino lumps at all masses in the mass range 
covered by our halo sample continuously grows from z — 1.^ until 
z = 1.4 (see the red, green, and blue curves in Fig.[8]representing 
z = 1.7, z = 1.6, and z = 1.4, respectively) reaching a maximum 
amplitude at 2; = 1.4. This is consistent with the interpretation 
given in Sec. |5.2| and exemplified in Fig.|7] between 2; = 1.8 and 
2; = 1.4 the neutrino mass grows in time making the neutrino po- 
tential wells progressively deeper, and at the same time the friction 
term acts as a real friction, slowing down neutrino particles and 
favoring the formation of bound structures, whose number density 
consistently grows in time. 

At 2; = 1.4, however, the scalar field (j) inverts its motion 
along its oscillations around the minimum of its effective potential, 
and consequently the neutrino mass starts decreasing, inducing the 
decay of the already formed neutrino potential wells. As a conse- 
quence of the inversion of the scalar field motion, also the friction 
term changes its sign, becoming an effective drag that accelerates 
the motion of neutrino particles, thereby increasing the neutrino ki- 
netic energy (see Fig. [TJ. This contributes to unbind a significant 
number of neutrino particles from previously bound structures, and 
determines the "evaporation" of a fraction of the neutrino lumps 
that had formed in the previous stage of neutrino mass growth. This 
phenomenon can be clearly identified in Fig. [8] where the purple 
and the orange curves, representing the neutrino halo mass func- 
tion at 2; = 1.2 and 2; = 1, respectively, show a progressively 
smaller number density of neutrino lumps at all masses with respect 
to the maximum amplitude of the neutrino mass function reached 
at2; = 1.4. 

Once again, we are witnessing an oscillatory behavior of the 
growth of large scale neutrino structures. Our results therefore show 
for the first time that neutrino structures forming in the context 
of Growing Neutrino cosmologies are not stable objects but rather 
experience subsequent phases of growth and decay, driven by the 
background evolution of the cosmon. We call this behavior "pulsa- 
tion" of the large scale neutrino structures. 

The detailed investigation of the "pulsation" of neutrino lumps 
and of their associated gravitational potential, combined with an 
estimate of the abundance of neutrino structures as a function of 
mass at different redshifts as provided by Fig. [8] will allow a deter- 
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mination of the evolution of the overall cosmological gravitational 
potential in Growing Neutrino cosmologies. 

5.5 Drag of Cold Dark Matter structures 

As a final step in our analysis, we now try to estimate the impact 
that the coupling between the scalar field and massive neutrinos 
can have on the evolution of CDM large scale structures. Since no 
coupling between DE and CDM is present in the specific scenario 
under investigation in the present work, the only way in which our 
model can alter the evolution of structure formation for the CDM 
density field is through a backreaction effect of the large neutrino 
structures that form below Znv onto the dynamics of CDM parti- 
cles. In particular, the formation of large neutrino structures at high 
redshifts can generate significant gravitational potential wells on 
large scales that might have a measurable impact on the properties 
of CDM large scale structures. 

In order to quantify this impact, we will make use of the three 
ACDM simulations described in Table [T] to directly compare the 
evolution of CDM structures in the presence of a coupling between 
DE and massive neutrinos to the standard ACDM evolution. 



5.5. 1 The CDM power spectrum 

First of all, we compute the CDM density power spectrum P{k) as 
a function of wavenumber k in all our simulations, and we compare 
the results obtained from simulations of the same scale with and 
without a coupling between DE and massive neutrinos. The results 
are shown in Fig.|9] where the ratio P(/c)/Pacdm (k) between the 
CDM power spectrum computed with and without neutrino cou- 
pling is shown for three different redshift values, namely z — 1.7, 
z = 1.3, and z = 1, for the smallest (left panel) and the biggest 
(right panel) box in our simulations set. 

As Fig.|9]shows, aiz = 1.7 there is basically no difference be- 
tween the two power spectra in any of the shown simulation pairs. 
This indicates that at this early time, the low neutrino mass and the 
small number of neutrino structures present in the universe are not 
sufficient to influence the dynamics of CDM structures, and that at 
this stage CDM is still acting as a seed for the development of the 
first neutrino overdensities. 

At z — 1.3, instead, an increase of CDM power of the or- 
der of ^ 1% at all scales in the simulations of Growing Neutrino 
cosmologies as compared to ACDM appears clearly in both plots. 
The largest scale simulations (ACDM-c and Cnu-c, right panel of 
Fig. [9} also show a weak scale dependence of the power enhance- 
ment, that appears to be larger at the largest scales of the simulation. 
This shows that dX z = 1.3 the gravitational potential induced by 
the formation of large neutrino structures has already a small but 
non-negligible impact on the underlying CDM distribution. The 
roles of these two cosmic component are therefore inverted: while 
at earlier times the CDM drives the motion of massive neutrino par- 
ticles, acting as a seed for the formation of the first neutrino lumps, 
now it is the fast growth of neutrino structures to mildly influence 
the evolution of the CDM distribution on very large length scales. 

The effect becomes more pronounced at z = 1, where the in- 
crease in CDM power can reach a level of 3 — 4% for wave numbers 
smaller than /c ^ 0.1 h/Mpc (left panel of Fig. [5]). At these scales, 
the ratio of power spectra shows significant irregularities. These are 
due to the statistical limitation determined by the box size, where 
only a few neutrino lumps are still present at z = 1 in the 40 (as 
well as in the 120 Mpc/h) simulation boxe (see Fig.|4]). This clearly 



induces an excess of power for some specific scales (determined 
by the specific spatial distribution of the lumps) that is reflected 
in the irregular shape of the power spectrum ratio in the left panel 
of Fig. [9] Our largest coupled neutrino simulation, however, with 
a size of 320 Mpc/h still contains a significant number of neutrino 
lumps at z = 1 (see Fig.|6}, thereby providing a statistically more 
reliable estimate of the CDM power spectrum. For this simulation 
(right panel of Fig. [9]), in fact, we see that no significant irregularity 
is present in the ratio of the power spectra with and without grow- 
ing neutrinos, and we can more easily quantify the backreaction 
of neutrino structures on CDM. As the plot shows, there is a clear 
scale dependence of the power enhancement, that reaches values 
of the order of ^ 5% at the largest scales of the simulation, while 
never exceeding 2% for scales smaller than /c ^ 0.1 h/Mpc. 

The backreaction effect of the formation of neutrino halos on 
the CDM large scale structures is therefore a potentially observable 
feature of the model, especially at very large scales, and could con- 
tribute to significantly constrain the Growing Neutrino scenario. 

At this point we also mention another characteristic feature of 
structure formation in Growing Neutrino quintessence. Measure- 
ments of the power spectrum by distributions of galaxies presum- 
ably give access to the CDM power spectrum, if we assume that 
baryons trace CDM (up to bias). On the other hand, for a measure- 
ment of the power spectrum by weak lensing the propagation of 
photons is influenced by the total gravitational potential, including 
the one induced by neutrino lumps. The power spectrum "seen" by 
photons is a spectrum for fluctuations of combined CDM, baryons 
and neutrinos. It may differ from the CDM spectrum probed by 
galaxy distributions, leading to inconsistencies between the ob- 
served galaxy and lensing power spectra if a standard ACDM cos- 
mology is assumed as the true underlying model. 

5.5.2 The CDM bulk flow 

As a last step in our analysis we consider the impact that the growth 
of large scale neutrino structures can have on CDM peculiar veloc- 
ities. We compute the average CDM peculiar velocity for all our 
simulation snapshots and we compare its time evolution for the 
coupled neutrino simulations with the corresponding ACDM runs. 
The results of this analysis are shown in the three plots of Fig. [T0| 
for the three different box sizes of our simulations, namely 40, 120, 
and 320 Mpc/h in the left, center, and right plots, respectively. 

All the three plots show a clear increase of the average CDM 
peculiar velocity at redshifts where the large scale neutrino struc- 
tures start playing a role in the coupled neutrino cosmologies. As a 
consequence of the sudden development of very large scale gravita- 
tional potential wells determined by the growth of neutrino lumps, 
CDM particles receive an additional acceleration with respect to 
the standard ACDM scenario, and progressively acquire larger ve- 
locities. This effect, as shown in Fig.[T0] determines an increase of 
the average CDM velocity of the order of 5 — 10%, and shows a 
clear scale dependence, with the weakest impact (^ 5%) for the 40 
Mpc/h simulation, where only a few mid- size (5 — 10 Mpc/h) neu- 
trino lumps are present at the redshifts under consideration, and the 
largest impact (^ 10%) for the 320 Mpc/h simulation, where in- 
stead a large number of intermediate- to-large size (10 — 50 Mpc/h) 
neutrino lumps can source the additional gravitational potential re- 
sponsible for the velocity excess. 

We have therefore shown that an anomalously large CDM ve- 
locity field might indicate the presence of an unexpected large-scale 
gravitational potential and that such potential could be sourced in 
our model by large scale neutrino structures. Growing neutrino 
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cosmolgies might therefore provide a natural explanation to the 
claimed detected large CDM bulk flow on scales larger than ^100 
Mpc/h HWatkins et aL|2009| ). 

Our preliminary analysis involves only the average peculiar 
velocity of all the CDM particles in the simulations. For a direct 
comparison with observational data, however, a determination of 
the velocity field for visible tracers of the CDM distribution through 
the construction of a full mock halo catalog would be required. We 
defer this analysis to a follow up paper. 

It is interesting to notice that the growth of the neutrino grav- 
itational potential at z > 1.3 is sufficient to induce a 10% velocity 
excess with respect to ACDM on scales larger than 300 Mpc/h, 
despite the subsequent decay of the neutrino gravitational potential 
at 1.0 < z < 1.3. This clearly shows that a detected anomalous 
CDM bulk flow in the low-redshift universe might not imply the si- 
multaneous existence of anomalously large-scale gravitational po- 
tentials, but could be explained, as for the case of the Growing Neu- 
trino models investigated in this work, by an epoch of subsequent 
growth and decay of such potentials at high redshifts (z > 1) whose 
effects remain imprinted in the CDM velocity distribution. Indeed, 
if only gravity provides for the attraction of particles, the bulk flow 
and the density contrast are tightly correlated. Enhanced bulk flow 
combined with only a mild increase in the density contrast would 
require additional attraction beyond gravity ( ,Ayaita et al.,2009j , as 
realised in our model. 



As pointed out above, such a peculiar time evolution (that we 
defined as a "pulsation") of the large-scale gravitational potential 
might imprint distinctive features on the large-scale CMB angular 
power spectrum through the ISW effect, and a detailed investiga- 
tion of these possible features, which goes beyond the scope of this 
paper, will be a necessary test in order to assess the viability of the 
model. 



6 CONCLUSIONS 

In the present paper, we have described the results of the first N- 
body simulations of structure formation in the context of Grow- 
ing Neutrino cosmologies. Our numerical treatment allowed for 
the first time to follow the fully nonlinear evolution of a statisti- 
cally meaningful sample of neutrino lumps forming in a specific 
Growing Neutrino model up to 2; = 1. This limiting redshift is 
due to the rapid growth of neutrino velocities into the relativis- 
tic regime, which does not fulfill the condition of small velocities 
that defines the newtonian limit for which standard N-body algo- 
rithms are implemented. Substantial modifications of the presently 
available numerical codes will be required in order to include in 
the treatment relativistic corrections, so to extend the reliability of 
numerical simulations down to 2; = 0. Nevertheless, our results 
significantly improve previous analyses by following the nonlinear 
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evolution of neutrino lumps to lower redshifts and well beyond the 
onset of virialisation. Furthermore, the use of N-body simulations 
allows to include the effects of directionality in the neutrino drag 
force and of merging and interaction of a realistic distribution of 
neutrino lumps. 

Our analysis has highlighted for the first time the phenomenon 
of pulsation in the growth of the density contrast in the neutrino 
fluid. In all realistic models of Growing Neutrino quintessence an 
oscillatory change of the neutrino mass has been found to set in 
once the interaction with the neutrino fluid becomes important for 
the time evolution of the cosmon and therefore for dark energy. An 
oscillatory neutrino mass results in a periodic change between fric- 
tion and drag for the neutrino velocities. For periods of increasing 
neutrino mass an effective friction term accounts for the fact that for 
constant momentum the velocity is reduced due to the increase in 
mass. In contrast, for periods of decreasing neutrino mass the oppo- 
site "drag term" enhances the neutrino velocities. Furthermore, the 
strength of the attractive force between neutrinos is proportional to 
the squared neutrino mass and shows oscillatory behavior for oscil- 
latory mass. The basic ingredients for the pulsation, i.e. the subse- 
quent increase and decrease in the neutrino velocities and density 
contrast that we observe in our N-body simulations, seem to be 
rather robust. We therefore expect that the phenomenon of pulsa- 
tion will be found also in more accurate treatments of the nonlinear 
growth of density fluctuations, which may overcome the limitations 
of the present numerical investigation. 

A pulsation in the growth of structure does not occur in the 
standard ACDM scenario. This phenomenon could therefore lead 
to interesting tests of the Growing Neutrino model. In particular, 
the CDM bulk flow may be enhanced despite very moderate modi- 
fications of the density contrast as compared to the ACDM model. 
Also the tension between the size of observed correlations of the 
ISW-effect in CMB anisotropics with observed structures on one 
side, and the observed bulk flow on the other side, could possibly 
be reduced. Quantitative statements will need, however, a treatment 
of the nonlinear structure growth that goes beyond the limitations 
of the N-body simulations presented in this work. 

A second important observation of the present investigation 
concerns the overall strength of the neutrino induced gravitational 
potential. We defer the quantitative analysis of the fluctuations in 
the gravitational potential to a separate publication. Nevertheless, 
already at the present stage it becomes apparent that the gravi- 
tational potential does never become very large for the redshifts 
z > 1 and the size of the boxes investigated in this paper. This can 
be inferred from the very modest increase in the CDM power spec- 
trum shown in Fig|9] This qualitative finding of limitations in the 
growth of the fluctuations in the gravitational potential confirms the 
arguments advocated by Pettorino et al. (2010). It clearly demon- 
strates that an extrapolation of the linear growth of neutrino density 
fluctuations and the neutrino induced gravitational potential (see 
e.g. [Franca et al.|2009| ) leads to huge overestimates for these quan- 
tities. 

Our third finding concerns the large neutrino velocities that are 
induced by the strong mutual attraction of lumps due to the cosmon 
force. The observation in our simulations that neutrino velocities 
become comparable to the velocity of light, limits the range of va- 
lidity of our treatment. On the other hand, the effective pressure for 
the neutrino fluid that is induced by these velocities counteracts the 
strong attraction mediated by the cosmon force. This pressure will 
limit the formation of very large and deep gravitational potential 
wells induced by neutrino lumps. In turn, the quantitative under- 
standing of the evolution of fluctuations in the gravitational poten- 



tial is crucial for the ISW-effect in CMB anisotropics. A quantita- 
tive computation of the imprint of Growing Neutrino quintessence 
for the CMB will have to wait for a quantitative computation of the 
cosmological neutrino induced gravitational potential. 

Presumably, this will require an extension of our methods to 
include effects from relativistic velocities. We do not know yet to 
what extent such effects will reduce the strength of the pulsation in 
the range z ^ 1. 

Finally, a fundamental new result of the present work is the de- 
termination of the statistical distribution of neutrino halos as a func- 
tion of redshift. For the first time we have attempted a computation 
of the neutrino halo mass function and of the quantitative effects 
on the CDM power spectrum and bulk flow. These estimates are 
still limited in redshift range to z ^ 1, and substantial uncertain- 
ties remain from the validity of the use of non-relativistic equations 
and from neglecting "backreaction" effects as well as locally vary- 
ing neutrino masses. Nevertheless, interesting effects which may 
lead to observational tests for the model become visible. These con- 
cern, in particular, an enhanced CDM bulk flow and the presence 
of very large scale neutrino lumps and associated gravitational po- 
tential wells. 



ACKNOWLEDGMENTS 

This work has been supported by the TRR33 Transregio Collabo- 
rative Research Network on the "Dark Universe", and by the DFG 
Cluster of Excellence "Origin and Structure of the Universe". Sup- 
port was given to V.P. by the Italian Space Agency through the 
ASI contracts Euclid-IC (1/031/10/0). All the numerical simula- 
tions have been performed on the Linux Cluster at the LRZ com- 
puting centre in Garching. 



APPENDIX 

In this paper we have applied the standard Newtonian equation of 
motion to all particles in the simulation. However, as we have seen, 
neutrinos reach soon relativistic velocities and after z 1 the as- 
sumption of non-relativistic motion is no longer acceptable. For 
completeness and in view of future applications we write down 
the equations of motion in the weak and slow-varying field limit 
but without limitations on velocities. We also include a non-linear 
scalar field potential. 

We assume the usual longitudinal metric for scalar perturba- 
tions 

ds'^ = -(1 + 2^)dt'^ + a^(l + 2^)dx'dxi (14) 
Defining the affine parameter r such that dr^ = —ds^ we can write 
dr = dt{l + ^)r (15) 

where for small ^, 

r= [l-^;^(l + 2($-^))]'/^ 
and where we define the peculiar velocity 



We need to derive the geodesic equations for a particle 

^+T^^^u-u^=0 (17) 
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where the four velocity is 



dr ^ r 'ar^ '''' 
Here we assume that the gravitational potentials are small, but the 
velocity can be relativistic. The geodesic equation holds for an 
uncoupled particle; below we extend it to coupled particles. As 
usual in the weak-field limit, we assume from now on that the time 
derivatives of ^, $ are negligible. In order to keep the equations 
linear in the potentials, we need to assume also ^ 1 (and 
similarly for Notice that d/dt is a total derivative and there- 

fore 

- ^ -[V +v ^ '-—^^ -) 

(19) 

where v —\v\ . 

From the geodesic equation with /j. = iwQ obtain to first order 
in ^: 

avi{l + si) + 'j'^avivv{l + S2) — 

-aHvi{l + £1) + v^^,i - (7^ + l)vi{vjVj^) + 

-f\^{vJ\/J^)-{^,^-V^^,^) (20) 



with 

= 2^(1 - 2717^) 
and where we used the usual special-relativistic factor 

Now as long as dark matter density is the dominating component, 
we can assume that the anisotropic stress induced by the relativistic 
neutrinos is negligible and therefore we put $ = — ^. Then we 
obtain 

Vi{l + Sl) + 'j'^ViVv{l + S2) = 

-HV^{1 + £1) + (27' + 1)V^{VJVJ^) - (1 + V^)^,^ (23) 

and from now on all spatial derivatives are with respect to physical 
distances ri = axi. The /i = equation is 

7^^;i;(l + S2) = -v^H{l + ei-4^) + {2-f^ -3){vjVj^) (24) 

(this can also be obtained up to terms higher order in ^ by multiply- 
ing Eq. ( [23| ) by v^) and we can use it to simplify Eq. ( [23| l obtaining 

7^v(l + si) = -ifv(l-2^)+47V(vV^)-(l + 2i;^7^)V^ 

(25) 

Now, to generalize to the coupled model, we can derive the 
coupled equation by conformally transforming the geodesic equa- 
tion. If we put 



_ 2/30 



we obtain 



and 



(26) 



(27) 



(28) 



Then, by assuming also that cp is taken at first order and (p is negli- 
gible we have 

v(l + £i) + j'^VVv{l + £2) = 

-^v(l + £1) - V(^ - /3ip) - v^V{^ + /3ip) + 
v(v • V(^ + /3(p)) + 27V(v • V^) 



where H = H{1 — /3^) . This reduces to the standard case 

V = -^v - V(^ - I3cp) (30) 
for small velocities. The /x = equation is 

7^^;i;(l + £2) = -^^;^(l + £i-4^) + (vV((27^-3)^ + /3(^)) 
which again can be used to simplify Eq. ([29|: 



(31) 



7^v(l + £i) = 

-Hv{l - 2^) + 47V(v • V^) - V(^ - Pep) - 27^^;V^ 

(32) 

This can be further generalized by including the higher-order terms 
in the perturbed field ip: 

7^v(l + £1) = 

-^v(l - 2^) + 47V(v • V^) - V(^ - /3ip) - 



27 + 2/3^V(p 
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